function emp_dist=emp_fit(shape,scale,n,mean_cycle_time)

% n is the max number of cycles
N=1000; % number of random values used from original distribution
emp_dist=zeros(n,1);
sum=0;
for i=1:N

    ran=gamrnd(shape,scale);
    cycle=ceil(ran/mean_cycle_time);
    emp_dist(cycle,1)=emp_dist(cycle,1)+1;
    sum=sum+1;
end

for i=1:n
emp_dist(i,1)=emp_dist(i,1)/sum;
end  